close all;clear all;clc;
% Parameters
g = 2;
lambda = 0.2;
eta = 0.6;
gprime = 0.5;

% Functions
tauf = @(g, gprime, lambda, eta, A, k)...
    1/(g-gprime)*(log(A)-log((exp((g-gprime)*k*eta)-1)/((g-gprime)*eta)...
    + (1-k)*((g/lambda)*exp((g-lambda-gprime)*k*eta)...
    - (g-lambda)/lambda*exp((g-gprime)*k*eta))));

alphaf = @(g, gprime, eta, A, k, tau)...
    (A*exp(gprime*(tau+k*eta))...
    - exp(g*tau)*(exp(g*k*eta)-exp(gprime*k*eta))/((g-gprime)*eta))/(1-k);

Abar = @(g, gprime, lambda, eta)...
    (exp((g-gprime)*eta/(lambda*eta)*log(g/(g-lambda)))-1)/((g-gprime)*eta);

A0f = @(g, gprime, lambda, eta, ktilde)...
    ((exp((g-gprime)*ktilde*eta)-1)/((g-gprime)*eta)...
    + (1-ktilde)*((g/lambda)*exp((g-lambda-gprime)*ktilde*eta)...
    - (g-lambda)/lambda*exp((g-gprime)*ktilde*eta)))...
    /(exp((g-gprime)*ktilde*eta)-(g-gprime)*eta...
    *((exp((g-gprime)*ktilde*eta)-1)/((g-gprime)*eta)...
    + (1-ktilde)*((g/lambda)*exp((g-lambda-gprime)*ktilde*eta)...
    - (g-lambda)/lambda*exp((g-gprime)*ktilde*eta))));

hf = @(g, lambda, eta, k) lambda*exp(lambda*k*eta)/(exp(lambda*k*eta)-1);

ktildef = @(g, gprime, lambda, eta)...
    fsolve(@(k) (1-k)*g*eta + g/lambda - (g-lambda)/lambda*exp(lambda*k*eta)...
    - exp(-(g-lambda-gprime)*k*eta), 0);

knc = @(g, gprime, eta, A)...
    1/((g-gprime)*eta)*log(A*(g-gprime)*eta+1);

% Calculations
ktilde = ktildef(g, gprime, lambda, eta);
A0 = A0f(g, gprime, lambda, eta, ktilde);
step = 1e-3;
A = [1.05 1.15 A0 1.25];
k = NaN([length(A) round(1/step)+1]);

for i = 1:length(A)
    kend(i) = floor(1/((g-gprime)*eta)*log(A(i)*(g-gprime)*eta+1)*round(1/step))+1;
    k(i,1:kend(i)) = 0:step:1/((g-gprime)*eta)*log(A(i)*(g-gprime)*eta+1);
    for j = 1:length(k(i,:))
        if ~isnan(k(i,j))
            tau(i,j) = max(tauf(g, gprime, lambda, eta, A(i), k(i,j)),0);
            tauPlusKeta(i,j) = tau(i,j) + k(i,j)*eta;
        else
            tau(i,j) = NaN;
            tauPlusKeta(i,j) = NaN;
        end
    end
end

% Points for figures
for i = 1:length(A)
    for j = 1:length(k(i,:))
        if (j==1)|(j==length(k(i,:)))|(mod(j,100)==0)
            k_figure(i, j) = k(i, j);
            tauPlusKeta_figure(i, j) = tauPlusKeta(i, j);
        else
            k_figure(i, j) = NaN;
            tauPlusKeta_figure(i, j) = NaN;
        end
    end
end

% Figure
h = figure;
set(h, 'defaultLegendInterpreter','latex',...
    'defaultAxesTickLabelInterpreter','latex');
a = axes;

plot(k(1,:), tauPlusKeta(1, :),'-.k', 'linewidth', 2);
hold on;
plot(k(2,:), tauPlusKeta(2, :),'-.k', 'linewidth', 2);
plot(k(3,:), tauPlusKeta(3, :),'-k', 'linewidth', 2);
plot(k(4,:), tauPlusKeta(4, :),':k', 'linewidth', 2);

p1 = plot(k_figure(1,:), tauPlusKeta_figure(1,:),'-.k*',...
    'linewidth',2,'MarkerSize',8);
p2 = plot(k_figure(2,:), tauPlusKeta_figure(2,:),'-.ks',...
    'linewidth',2,'MarkerSize',8);
p3 = plot(k_figure(3,:), tauPlusKeta_figure(3,:),'-k',...
    'linewidth',2,'MarkerSize',8);
p4 = plot(k_figure(4,:), tauPlusKeta_figure(4,:),':kd',...
    'linewidth',2,'MarkerSize',8);

xl = xlim();
yl = ylim();
set(a, 'xlim', [0 1], 'ylim', [0 0.55], 'XTick',[0 0.5 ktilde 1],...
    'XTicklabel',{'$0$' '$0.5$' '$\tilde{k}$' '$1$'}, 'YTick',[0.1 0.3 0.5],...
    'YTicklabel',{'$0.1$' '$0.3$' '$0.5$'}, 'fontsize', 16);
plot([ktilde ktilde],[0 tauf(g,gprime,lambda,eta,A(end),ktilde)+ktilde*eta],...
    'k--', 'linewidth', 1, 'HandleVisibility','off');
aux0 = xl(1):0.02:xl(2);
auy0 = aux0*eta;
plot(aux0, auy0, 'k.','MarkerSize',4, 'HandleVisibility','off');

text(0.2, 0.1, 'slope $=\eta$', 'interpreter','latex', 'Fontsize', 16);

xlabel('$k$', 'fontsize', 20, 'interpreter','latex');
ylabel('$\tau^{*}+k\eta$', 'fontsize', 20, 'interpreter','latex');
lgd = legend([p1, p2, p3, p4], {'$A=1.05$', '$A=1.15$', '$A=A_0=1.20$', '$A=1.25$'},...
    'location', 'southeast', 'fontsize', 16);
lgd.ItemTokenSize(1) = 24;

%saveas(h, 'tauplusketa_A1.eps');


